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C*~) ' Abstract. To describe and simulate dynamic micromagnetic phenomena, we consider a 

coupled system of the nonlinear Landau-Lifshitz-Gilbert equation and the conservation of 
momentum equation. This coupling allows to include magnetostrictive effects into the sim- 
ulations. Existence of weak solutions has recently been shown in [12]. In our contribution, 
' we give an alternate proof which additionally provides an effective numerical integrator. 

The latter is based on lowest-order finite elements in space and a linear-implicit Euler 
time-stepping. Despite the nonlincarity, only two linear systems have to be solved per 
. timestcp. and the integrator fully decouples both equations. Finally, we prove uncondi- 

tional convergence — at least of a subsequence — towards, and hence existence of, a weak 
solution of the coupled system, as timestep size and spatial mesh-size tend to zero. Nu- 
1 merical experiments conclude the work and shed new light on the existence of blow-up in 

, micromagnetic simulations. 

rd 

<3 . 

1. Introduction 

^ ■ Throughout all technical areas, magnetic devices like sensors, recording heads, and magneto- 
resistive storage devices are quite popular and thus widely used. As their size decreases to 
a microscale, and the testing and development becomes more and more involved, the need 
for reliable and stable simulation tools as well as for a thorough theoretical understanding 
rises. In terms of mathematical physics, micromagnetic phenomena are modeled best by the 
. Landau-Lifshitz-Gilbert equation (LLG), see (1) below. This nonlinear partial differential 
equation describes the behaviour of the magnetization of some ferromagnetic body under the 
influence of a so-called effective field. The mathematical challenges as well as its applicability 
to a wide range of real world problems makes LLG an interesting problem for mathematicians 
and physicists, but also for scientists from related fields like engineers and developers from 
high-tech industry. 

In our contribution, we present and analyze a computationally attractive integrator to 
solve LLG numerically. Additionally, our analysis provides a constructive existence proof for 
weak solutions of the coupled system for LLG with magnetostriction and thus particularly 
includes the results of [12]. For the pure LLG equation, existence and non- uniqueness of 
weak solutions of LLG goes back to [3, 26] for a simplified effective field. For a review of the 
analysis of LLG, we refer to [13, 15, 22] or the monographs [20, 23] and the references therein. 
As far as the numerical analysis is concerned, mathematically reliable and convergent LLG 
integrators are found in [1, 2, 4, 5, 8, 9, 11, 16, 17, 19, 25]. Of utter interest are unconditionally 
convergent integrators which do not impose a coupling of spatial mesh-size h and time-step 
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size k to ensure stability of the numerical integrator. Those integrators are split into two 
groups; first, midpoint-scheme based integrators [4, 25] which rely on the seminal work [9] 
of Bartels & Prohl; second, projection-based first-order integrators [2, 5, 11, 16, 17, 19] 
which build on the work [1] of Alouges. All of the above integrators have in common that 
they allow for constructive existence proofs of the corresponding problems. The idea, which 
is also exploited in the current work, is to show boundedness of the computable discrete 
solutions. Then, a compactness argument concludes the existence of weakly convergent 
subsequences. Finally, those weak limits are identified as weak solutions of the coupled 
system. 

In our contribution, we extend the analysis of the aforementioned works for the projection 
based integrators and show that these ideas can also be transferred to the coupled system 
of LLG with magnetostriction. Even though the structure of the main proof in this work is 
similar to the one from e.g. [1, 5, 11] (and also the other works mentioned), we stress that 
the individual ingredients are much harder to gain than for the coupling with stationary 
equations [11] or with the instationary Maxwell system [5]. First, unlike the Maxwell-LLG 
system, the coupling to the conservation of momentum equation involves a nonlinear coupling 
operator. Second, the field contribution which accounts for magnetostriction, does not only 
depend on the displacement, but also on its spatial derivative. For these two reasons, the 
analysis requires new mathematical tools and therefore complicates the convergence proof. 
In contrast to the existing literature, we thus extend the approach from [1] and analyze a 
nonlinear coupling of LLG to a second time-dependent PDE. The contributions of this work 
as well as the advances over the state of the art can be summarized as follows: 

• We include the magnetostrictive field contribution into the numerical analysis as well 
as into the algorithm to account for elastic effects on a microscale. This allows to 
conduct more precise simulations in certain applications. 

• We give a new and constructive proof for the existence of weak solutions for LLG 
with magnetostriction by providing a numerical integrator which is mathematically 
guaranteed to converge to a weak solution of the coupled system of LLG with mag- 
netostriction. 

• The proposed integrator is unconditionally convergent towards a weak solution as 
soon as timestep-size k and spatial mesh-size h tend to zero, independently of each 
other. For midpoint-scheme based integrators [4, 9, 25], unconditional convergence 
is theoretically proved. However, the solution of the nonlinear system requires either 
a heuristical solver or an appropriate fixed-point iteration. The latter is used in [4, 
9, 25], but the resulting explicit scheme again involves a coupling of h and k to 
guarantee convergence and avoid instabilities. 

• The proposed algorithm is extremly attractive from a computational point of view: 
Since the two equations are fully decoupled, the implementation is easy and only 
two linear systems have to be solved per timestep. Contrary, midpoint- scheme based 
integrators [4, 9, 25] have to solve one large nonlinear system per timestep, and the 
decoupling has not been thoroughly analyzed [4]. 

• We provide a numerical comparison between the extended Alouges-type integrator 
and an integrator based on the midpoint scheme [25]. Empirically, the numerical re- 
sults show that our integrator works with larger timesteps than the midpoint scheme. 
In particular, in our setting the computations based on the proposed integrator turn 
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out to be much faster. This gives empirical evidence that the fixed-point iteration 
for the midpoint scheme indeed poses a computational bottleneck. 

Outline. The remainder of this paper is organized as follows: In Section 2, we state the 
mathematical model for the coupled system of LLG with magnetostriction and recall the 
notion of a weak solution (Definition 1). In Section 3, we collect some notation and prelimi- 
naries, as well as the definition of the discrete ansatz spaces. In Section 4, we write down our 
numerical integrator in Algorithm 2, and Section 5 is devoted to our main convergence result 
(Theorem 4) and its proof. Finally, numerical examples concludes the work in Section 6. 

2. Model Problem 

The evolution of the magnetization of a ferromagnetic body Q during some time interval 
(0, T) is mathematically modeled by the Landau-Lifshitz-Gilbert equation (LLG) which in 
dimensionless form reads 

m t — am x m t = — m x H cff in Q T := (0, T) x Q (la) 

Here, m : §> 2 : = {x G M 3 : |x| = l} denotes the sought magnetization, and H e g- is the 

so-called effective field that consists of several energy contributions each of which models a 
certain micro magnetic effect. More precisely, in our work, the effective field consists of the 
exchange contribution Am, the magnetostrictive component h m , and all other stationary 
and lower-order effects are collected in some general field contribution 7r(m), i.e. 

H eff = C , e Am + h m -7r(m). (lb) 

The general energy contribution 7r(m), is only assumed to fulfill a certain set of properties, 
see (20)-(21), and we emphasize that this particularly includes the case H e g = C e Am + 
h m + Ca,niD&(m) + P(m) — f. Here, $(m) denotes the crystalline anisotropy density and 
f is a given applied field. The contribution P(m) stands for the nonlocal strayfield. The 
constants C e , C mi \ > denote the exchange and anisotropy constant, respectively. To keep 
the presentation simple, we did not include the additional coupling to the full Maxwell 
equations as in [5]. We stress, however, that this extension is straightforward and, with the 
combined techniques from this work and [5], one could consider a coupled system of full 
Maxwell LLG with the conservation of momentum equation to account for magnetostrictive 
effects. The combined results from [5] and the current work would then directly transfer to 
the coupled case and we would still derive unconditional convergence. Moreover, with the 
techniques from [11], it is straightforward to rigorously include a numerical approximation 
7T/ l (-) of 7r(-) into the convergence analysis. 

For a bounded Lipschitz domain Q C IR 3 and a time interval (0,T), we now aim to 
solve (la) on Qt supplemented by the initial and boundary conditions 

m(0) = m° G H 1 ^; § 2 ) and d n m = on (0, T) x dQ. (lc) 

The constraint |m°| = 1 almost everywhere in Q models the fact that we consider a constant 
temperature below the Curie point. Multiplication of (la) with m yields <9 t |m| 2 = 2m m t = 
0. Hence, the modulus |m| = 1 is preserved in time. By imposing the modulus constraint 
on m°, this guarantees |m(t)| = 1 for almost all times t G (0,T). 
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For modeling the magnetostrictive component, we follow the approach of Visintin [26]. 
Here, the magnetostrictive field reads 

3 

h m : ft r -)■ M 3 , (h m ) ? = (h m (u, m))^ := ^ \™ pq a tj (m) p , (2) 

where (-) p denotes the p-th component of a vector field. We implicitly assume linear depen- 
dence of the stress tensor cr = {cr^} on the elastic part of the total strain e e = {etj} which 
is the converse form of Hook's law, i.e. 

3 

cr := A e £ e (u, m) : ft T — ► M 3x3 , a %] = £ A^, (3) 

p,q=l 

£ e (u, m) := e(u) - e m (m) : [0, T] x ft — > M 3x3 , (4) 

where u : ft^ — y M 3 denotes the displacement vector field. The fota/ strain is defined by the 
symmetric part of the gradient of u, i.e. 



2 \<9xj _ 
and the magnetic part of the total strain by 

3 

e m (m) := A " mm' : ft T — ► M 3x3 , e £(m) = £ A™ ,(m) p (m),. (6) 

p,q=l 

In addition, we assume both material tensors A G {A e , A m } to be symmetric {\ij pq = Aji pg = 
\jqp — \qij) an d positive definite 

3 3 
^ ^ijpqiijipq > A* ^ £?• (7) 

ij,p,q=i *J=i 

with bounded entries, i.e. there exists some A with A? p? , A™ pg < A for any i,j,p, q = 1,2, 3. 
The stress tensor cr and the displacement field u (where we assume no external forces) are 
finally coupled via the conservation of momentum equation 

QVLu — V ■ cr = in fty. (8a) 

Here, we assume the mass density g > to be constant and independent of the deformation. 
Equation (8a) is additionally supplemented by the initial and boundary conditions 

u(0) = u° in ft, u t (0) = u° in ft, and u = on 9ft. (8b) 

Altogether, we thus aim to solve the coupled problem 

m t -amxm« = -mxKU 
gu tt - V ■ cr = 0, 

subject to the stated initial and boundary conditions. 
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Using the above boundary conditions, Hook's relation (3), the definition of the total strain 
tensor (4), and the symmetry of the tensors A e and A m , we obtain the following variational 
formulation of (8a) 

(gu tt (t),cp) + (\ e e(u)(t),e(cp)) = (AV"(m)(t), e(<p)) for all <p G Hj(fi). (10) 

Given these notations, we now define our notion of a weak solution for the coupled LLG- 
magnetostriction system (9), which is the same as in [12]. 

Definition 1. The tupel (m, u) is called a weak solution of LLG with magnetostriction, if 
for all T > 0, 

(i) m G H 1 (f2j') with |m| = 1 almost everywhere in Qt o,nd u G H 1 (Q,t); 

(ii) for all 4> G C°°(tt T ) and C e ([0, T); C°°(ft)) , we fcaue 

/ (m t ,4>)-a ((m x nii)» = -Ccxch / ((Vm x m), V0) 

+ [ ((h m xm),0)- / ((ir(m)xm),0) (11) 

-ff/ f (A e e(u), £ (C)) = / (A e s(m),s(C))+ /(u°,C(0,-)) (12) 

(iii) t/iere ZioWs m(0, •) = m° and u(0, ■) = u in the sense of traces; 

(iv) for almost allt' G (0,T), we nave bounded energy 

l|Vm(OH^ ( n) + llmtH^j + ||Vu(OllL«(n) + K(OllL«(n) < d, (13) 
where C\ is independent oft and depends only on |Q|,mo,Uo, and uq. 

3. Preliminaries 

For time discretization, we impose a uniform partition = to<ti<...<tN = Tof the 
time interval [0, T]. The timestep size is denoted by k = kj := tj+\ — tj for j = 0, . . . , N — 1. 
For each (discrete) function which is continuous in time, ipi = <p(tj) denotes the evaluation 
at time tj. For the time derivatives in the conservation of momentum equation (10), we use 
difference quotients of first and second order which are denoted by 

d t z i = , d tZi = = . (14) 

For spatial discretization, let T~h be a quasi-uniform, regular triangulation of the polyhedral 
bounded Lipshitz domain Q C 1R 3 into tetrahedra. The spatial mesh-size is denoted by h. By 
S 1 (7h), we denote the lowest-order Courant FEM space of globally continuous and piecewise 
affine functions from Q to M 3 , i.e. 

S\T h ) := {<f> h G C(Q;R 3 ) : <f> h \ K G V 1 (K) for all K G %}■ 

By 1^ : C(Q; M 3 ) — > S l (Th), we denote the nodal interpolation operator onto this space. The 
set of nodes of the triangulation Th is denoted by Afh- 

For the discretization of the magnetization m in the LLG equation (la), we define the set 
of admissible discrete magnetizations by 

M h := {<j> h G S\%) : \<j> h {z)\ = 1 for all z G M h }. 
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Furthermore, for 4> h G Mh, let 

:= tyft e S\T h ) : i^( z ) ■ ^(z) = for all z G A4} 

be the discrete tangent space associated with tf> h . Here, x • y stands for the usual Euclidean 
scalar product of x, y G M 3 which is sometimes also denoted by (•, •) to improve readability. 
The L 2 {Vt) scalar product is denoted by (•, •) throughout. Due to the modulus constraint 
m t • m = 0, and thus |m(£)| = 1 almost everywhere in Qt, we discretize the time variable 
v(ij) := m t (i,) in the discrete tangent space of m^. 

To discretize the equation of magnetoelasticity (10), we employ Sq (%) '■= 5 1 (7/ l ) C\Hq (fi). 
In addition, let m° G M.h an d u°, u° G SliTh) be suitable approximations of the initial data 
obtained e.g. by projection. Further requirements on those initial data are specified below 
in Theorem 4. Finally, we define d t u° as u°. Througout this work, we write A < B if there 
holds A < CB for some h and k independent constant C > 0. 

4. Algorithm 

For discretization of the LLG equation, we follow the approach of Alouges [1] which has 
been generalized in [2] and Goldenits et al. [11, 18, 19]. The main idea is to introduce 
a new free variable v w m t and to interpret LLG as a linear equation in v. This ansatz 
exploits the formulation 

am t + m x m t = H cff — (m ■ H cff )m, (15) 

which is equivalent to (la) under the constraint |m| = 1 almost everywhere, see e.g. [18, 
Lemma 1.2.1]. The conservation of momentum equation is discretized in space by a standard 
FEM approach, and by finite differences in time. This approximation is in analogy to the 
work of Banas and Slodicka [6], where the focus is on FEM discretizations for strong 
solutions of (1) with (8). In addition and for computational ease, the two equations can be 
decoupled. We propose and analyze the following algorithm: 

Algorithm 2. Input: Initial data mj and parameter < 9 < I, a > 0. For t = 
0, . . . , N — 1 iterate: 

(i) Compute unique solution \ l h G JC m e such that for all tp h G JC m e , we have 

(ii) Define m e h +1 G M h nodewise by m£ +1 (z) = gr|^f^yj for all z G Af h . 

(iii) Compute unique solution u^ +1 G S^Th) such that for all if? h G S^Th), we have 

g(4n{ +1 ^ h ) + (A^u^ 1 ), £(</>,)) = (A e s m K +1 ), £(</>,)). (17) 
In the above algorithm, the discrete magnetostrictive contribution is given by 

3 

[h m (ui,m£)] 3 := Yl X m AM) P > with <r h = * e {e(u e h )-e m (mi)). 

i,3,P=l 

Exploiting that we solve each of the two equations separately, we can immediately state 
well-posedness of Algorithm 2. 
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Lemma 3. Algorithm 2 is well defined, i.e. it admits unique discrete solutions (v£, m^ +1 , u^ +1 ) 
in each step £ = 0, . . . ,N — 1 of the iteration. Moreover, we have ||m^||i,°°(n) = 1 f or oil 

e = i,...,N. 

Proof. We first show solvability of (16). We define the bilinear form 

a{(; •) : JC mi x K mi M, a{(<j>, if) := a(<p, tp) + 9C e k(V(j>, Vc^) + (K x cp), <p) 

and the linear functional L\(ip) := C e (VmJ, Vy>) + (h m (u£, m£), <£>) — (7r(m|),<^). Then, 
(16) is equivalent to a{(v e h ,ip h ) = L\(ip h ) for all ip h G /C m *. Note that a[(-,-) is positive 
definite for a > 0, i.e. a e ((p,(p) > || V 5 II l 2 (s^) • Thus, by exploiting finite dimension, we see 
that there exists a unique v e h G /C m ^ which solves (16). Due to pointwise orthogonality of m e h 
and v£, and the Pythagoras theorem, we get |m^(z) + kv e h (z)\ 2 = |m^(z)| 2 + k\v e h (z)\ 2 > 1 
and thus even step (ii) of the above algorithm is well-defined. The bound ||m^||L°°(n) = 1 can 
be seen by the normalization at the grid points in combination with barycentric coordinates 
and the convexity of each tetrahedron. 

For the second equation (17), we consider the bilinear form 

aa(v) : S]{T h ) x S]{T h ) ->R, a 2 (C, V) := t|(C,</0 + (A e e(C),^W) 

and the linear functional L e 2 (if>) = (A e £ m (m^ +1 ), £ (*/>)) + §(dtu£, i/>) + p(u^, According 
to (7) and Korn's inequality [10, Thm. 11.2.16], it holds that a 2 (ip,if>) > ||V>IIhw With 
this notation, (17) is equivalent to a2(u^ +1 , i/j h ) = L^fyh) for all functions ip h G S^Th), and 
hence, admits a unique solution u^ +1 G <So(Th) in each step of the loop. □ 

5. Main Theorem 

In this section, we aim to show that the preceding algorithm indeed converges towards the 
correct limit. Before we start with the actual analysis, we collect some general assumptions 
and some more notation. Throughout, we assume that the spatial meshes Th are uniformly 
shape regular and satisfy the angle condition 

/ VCi • VCi < for all hat functions Q, Q G S l (T h ) with i ^ j. (18) 
Jq 

This somewhat technical condition is a crucial ingredient of the convergence proof, since it 
yields the discrete energy decay 

llVm^H^ < ||VK + kvi)\\l m . (19) 

The estimate (19) is a direct consequence of the inequality || VI/i(m/|m|) H^arm < || VX^mH^^) 
which was proved by B artels in [7]. We like to emphasize that the angle condition is al- 
ways fulfilled for tetrahedral meshes with dihedral angles that are smaller than tt/2, cf. [7], 
and easily preserved for uniform mesh- refinement. For x G Q and t G f^,^+i) and for 
7^ G {m e h ,v e h ,u e h }, we define the time approximations 

y*{t,x) :=^7r(x) + ^^(x), 
7,- fc (t,x):= 7 I(x), 7, + fc (t,x):= 7 r(x). 
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Note that 77^ can also be written as 7/ l fc(£, x) = 7^(x) + (t — t£)dt7^ +1 (x). In addition, for 
t G we define 

uwfcftx) := d t u£(x) + (t - t,)d 2 < +1 (x), u^(t,x) := d t u£(x), u+ (t,x) := d t < +1 (x). 

The next statement is the main theorem of this work and particularly includes the main 
result from [12]. 

Theorem 4. (a) Let 9 G (1/2,1] and suppose that the meshes Th are uniformly shape 
regular and satisfy the angle condition (18). Moreover, let the general energy contribution tt 
be uniformly bounded in L 2 (Qt), 

IKWIIl 2 ^) — f or a M l n l e L 2 (f2 T ) zi/ii/i n < 1 a.e. in (20) 

u>rf/i an n-independent constant C n > anc? assume weak convergence of the initial data, 
i.e. — ^ m°,u° — ^ u° in H 1 (f2) ; as well as u° — u° in L 2 (f2) as h — > 0. Under 
these assumptions, we have strong L 2 (f2^)- convergence of m^ k towards some function m G 

(b) Suppose that, in addition to the above assumptions, we have 

7r(m^ fc ) — 1 7r(m) weakly subconvergent in L 2 (f2r)- (21) 

Then, the computed FE solutions (m^, u^) are weakly subconvergent in H 1 ^?-) x H^f^) 
towards some functions (m, u), and those weak limits (m, u) are a weak solution of LLG 
with magnetostriction. In particular, weak solutions exist and each weak accumulation point 
of (mfcfc, U/jfc) is a weak solution in the sense of Definition 1. 

The proof will be done in roughly three steps. 

(i) Boundedness of the discrete quantities and energies. 

(ii) Existence of weakly convergent subsequences. 

(iii) Identification of the limits with weak solutions of LLG with magnetostriction. 

As mentioned, we first show the desired boundedness and start with some preliminary lem- 
mata. 

Lemma 5. The discrete magneto strictive component can be estimated by the total strain of 
the discrete displacement, i.e. 

||h m «,mi)|| 2 L2(n) < C 2 ||6«)|| 2 2(n) + C 3 , (22) 
for some constants C2, C3 > that depend only on A. 
Proof. By definition of the magnetostrictive part, we immediately get 

HilmtUjj, Ulh)\\lP(p) ~ A ll m /illL°°(n)ll <7 llL 2 (Q) ~ ll CT llL 2 (Q) 

due to the normalization step. From the definition of the discrete stress tensor and the 
boundedness of the material tensors, we additionally get 

\\<r h \\h m < M\e(ui)\\l (n) + ||£ m K)|| 2 2(n) < Mvfowiw + c. 

This yields the assertion. □ 
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Lemma 6. For j = 1, . . . , N, there holds 

l|Vm{|£ 2(n) + (6- l -)k* Y}\ VviHL(o) + *E llvill^n) 

.=o ^ (23) 

< Gi(||VmS||ia (n) + kJ2 \\<ni)\\l 2{n) + C 5 ), 

1=0 

for constants C4, C5 > t/iat depend only on C„ as well as C2 and C3 from the previous 
lemma. 

Proof. In (16), we use the special test function ip h = v e h G /C m £ and get 

a(vj,v£) + ( K x yj),vp = -C e {V(m{ + ekvi),Vvi) + (h m (u£, m£), v£) - (7r(m£)X) : 

=0 

whence 

«l|v£|&<n, + Ce0*||Vvi||k (n) = -C e (Vml, Vv£) + (h m (4,ml),va - (7rK),vl). 
Next, we use the fact that || Vm^y 2 ^^ < ||V(m£ + fcv^)!]^^)) see (19), to obtain 

^l|VmJ +1 ||i 2(n) < i||VmX (n) + fc(Vmi, Vvi) + y || Vv£|| 2 2(n) 

< ^l|Vmi||^ (0) - (9 - 1/2)^1^^11^^ (24) 

- ttKIIl^q) + ^-(h m «,m^), V y - -^(7r(m£), v£) . 
We sum over the time intervals from to j — 1, which yields for any v > 

1 afc j ~ 1 

2ll Vm hllL 2( n) + ^rEl|v[||^ { n) 

e £=0 



< ^HVmX (n) + ^E[(h m (u£,m£),v£) - (7r(mi),vi)] - (61 - l/2)fc 2 £ II Vv£||2, a(n) 



£=0 £=0 



< ^<\\h { n) + ^E (Hhm(ul,ml)||S, 2(n) + ||ir«)||k<n,) 

e £=0 

+ ^ E II vj 112,(0) - (5 - l/2)fc 2 E II Vvil|^ 2(n) =: RHS. 

e £=0 £=0 

With Lemma 5 and the uniform boundedness (20) of 7r(-), we further obtain 
RHS < -||V m °||^ (n) + — 1 E (lk(ul)|| 2 L2( n) + C 3 ) + C w 



2 n /.ni^u; 4C 

+ ^ E 11 v ^Il^) - {0 - 1/2)^ E 11 v<\\t m - 

e 1=0 e=o 
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In total we thus derive 

j'-i , j'-i 



^llVmill^o, + £(a - i/) £ living + ((? - i)* 2 £ II Vv * 



|2 

'/illL 2 (Q) 

< ^|VmX (n) + |^£ + C 3 ) + a 

Taking v < a thus yields the desired result. □ 

Given the last two lemmata, we now aim to show boundedness of the discrete quantities 
involved in equation (17), i.e. boundedness of the discrete displacement approximations. The 
following result has basically already been stated in [6, Lemma 3]. Since we require some 
important modifications, we include the proof here. 

Proposition 7. For any j — 1, . . . , N and | < 9 < 1, there holds 

j j 
\\dt<\\h { n) + W^h ~ ^ u Ml^) + H £ «)IIl^) + £ H £ ( u £) _ e « -1 )llL3(n) < C 6 

(25) 

for some h and k independent constant C§ > which depends only on X* and the constants 
C4 and C5 from Lemma 6. 

Proof. We use ij) h = u^ +1 — u e h as test function in (17) and sum up for t = 0, . . . , j — 1 to see 
,(d^ +1 - d 4 uU< +1 ) + (A e e« +1 ),s« +1 ) - e (u£)) = (A^(m[ +1 ),sK +1 ) - e(u{)). 
Recall that the Abel summation formula shows for any «;£l and j > 
j'-i i-i 



e=o i=o 
This in combination with the symmetry and the positive definiteness (7) of A 6 now yields 

i J 

ll<Wj£ 2(n) + £ IK< - dA^Whin) + \\ £ (<)\\h(u) + Yl H £ K) - £ ( u i" 1 )llL2 ( n ) 

1=1 

< C + C£ (AV"K), «) ~ «% 



i=i 



i=i 



for some generic constants C,C > 0. Note that the terms Hdtu^H^^) = ||u°||l2(q) and 
||£(u°)||l2(q) are uniformly bounded due to the assumed convergence of these initial data 
and are hidden in the constant C . 
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Next, we rewrite the sum on the right-hand side as 



1=1 

3-1 

= (A e s m (mi),eK)) - (AV»K), e (u°)) - £ (A e e m K +1 ) - AV"K),s(u£)) 

e=i 

= (A e s m (mi),eK)) - (AV"K),s«)) - kf^ (A e d,6 m K +1 ),6«)). 

£=1 

For any 77 > and due to boundedness of A 6 , we further get 

3 3 

+ E ll d X - dtuj- 1 !!^^) + ||s«)||2 2(Q) + ^ || e ( u £) - eK- 1 )!!^^ 



^=1 

< 1 + kJ2 l|d*e m K+ 1 )|£ a(n) + ||e(ui)|fc (n) + ^lk m K)||^ (c) 



=1 £=1 



+ *?lkWII^(n) + lk m M)llL» (n) + IKOIIlw 
For sufficiently small 77, this can be simplified to 

\\dtu{\\l Hn) + E \\<k< - d^wi^ + Mufowhw + E lk(ui) - sK- 1 )!!^ 
£=1 £=1 

< 1 + k( ]T \\d t e™(m{ +1 )\\l {n) + £ || £ (u()||^ (n) ). 



Here, we again used convergence of the initial data. Using the boundedness and symmetry 
of A ,n in combination with boundedness of ||m£||L°°(n), straightforward calculation shows 
||d t e m (n4)||2 2(n) < ||d t ml||2 2(n) , c f. [24]. In combination with 

1 i—\ 

11 , £ 1 1 2 II 1 1 2 << 1 1 £ 1 1 1 2 

llw^IlL 2 ^) = II 1 \\L 2 (n) — ll v ft Hl 2 (Q); 

cf. [1, 2, 11] resp. [18, Lemma 3.3.2], this results in 

3 3 
I|dt4||^ (n) + E ||d t u£ - dtuj- 1 !!^^) + ||e(u£) + E ||e(u£) - eCu^ 1 )!^ 

£=1 £=1 
j'-i i-i 

<i+KEii v "X ( n)+Eii £ ( u i)llLW- 



£=1 
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Next, we apply Lemma 6 to see 

3 3 

l|dt<||^ 2(n) + l|d*u{ - dtul'm^ + ||e(u£) + E ||e(ul) - ^u^ 1 )!)^ 

£=1 1=1 

< 1 + (II ^<\\h { n) + k E |k«)||2 2(n) + fc £ ||£«)IIlW 

fcO £=1 

i-i 



<i + *5>(u£)||k (n) . 

£=0 

Application of a discrete version of Gronwall's lemma finally yields the assertion. □ 

In order to show the desired H 1 (f2r)-convergence of u, we still need to show uniform 
boundedness of the L 2 (f2^)-part. This result is stated in the next corollary. 

Corollary 8. Due to the boundary conditions employed, the Poincare inequality in combi- 
nation with Korn's inequality shows 

HWbm + E K-uMw < c 7 (\\vui\\l m + \M< - K^Whm) 

i=\ e=i 

< C 7 (\\e(u{)\\l {n) + £ || £ (u| - uj" 1 )!!^) (26) 

i=i 

for any j = 1, . . . , N, where the constant C 7 > stems from Poincare 's inequality and thus 
only depends on the diameter ofQ. According to Proposition 7, the right-hand side of (26) 
is uniformly bounded. □ 

Proposition 7 now immediately yields boundedness of the discrete magnetizations. 

Corollary 9. For any j = 1, . . . ,N, there holds 

\\^<\\h(U) + kJ2 \\<\\b(a) + ( e ~ \> 2 E H Vv £Hl^) < C 8 (27) 

1=0 1=0 

for some constant C 8 > which depends only on C 4 , C 5 , and C 6 . 
Proof. From Lemma 6, we get 



i-i j'-i 



|Vmi|| 2 L2(n) + (6- hk 2 Ell^||^ ( n) + k E llv* 1 



2 J / Ji v v h\W(n) -i~ » / j II v feiiL 2 (Q) 

1=0 1=0 

j-i 

< C 4 (||Vm°||^ (n) + fc^||e(ul)||^ (n) + C 5 ) ; 

1=0 



(2t 



By utilizing Proposition 7, we see &X^=o || £ ( u DllL 2 (n) — I ^ | C 6 . The uniform boundedness 
of ||Vm° ||i, 2 (n) concludes the proof. □ 

Next, we deduce the existence of convergent subsequences. 
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Lemma 10. Let 1/2 < 9 < 1. Then, there exist functions (m, u, u) e H 1 (f2 T , § 2 ) xH 1 (r2 r ) x 
L 2 (f2y) such that 

m.hk — nx m H^fijO (29a) 

m hfc , m± -mm ^(H 1 ), (29b) 

m^,, ->mm L 2 (fi T ) (29c) 

m^fc, m^ fc — )■ m pointwise almost everywhere in (29d) 

Um u in H^fir) (29e) 

u^u^uml^H 1 ), (29f) 

Uhfc, uj^. -^um L 2 (fi T ) (29g) 

Uhfc, u^fc -^iim L 2 (fi T ). (29h) 

i/ere, i/ie convergence is to be understood for a subsequence of the corresponding sequences 
which is successively constructed, i.e. for arbitrary spatial mesh-size h — > and timestep-size 
k — y 0, there exist subindices h n ,k n , for which the above convergence properties are satisfied 
simultaneously. In addition, there exists some v G L 2 (f2 r ) with 

v u -vmL 2 ((],) (30) 

again for the same subsequence as above. 

Proof. From Proposition 7, Corollary 8, and Corollary 9, we immediately get boundedness 
of all of those sequences. A compactness argument thus allows to successively extract con- 
vergent subsequences. Therefore, it only remains to show, that the corresponding limits 
coincide, i.e. 

\imj hk = lim7^ = lim7+ with j hk e {m hk , u hk , u hk }. 

To see this, we make use of the boundedness of two subsequent solutions. From Proposition 7, 
Corollary 8, and Corollary 9, we see that 

i-i j'-i 

\\ m h m AllL 2 (tt) + 2-^i 'I h U /'JlL 2 (Q) 

1=0 i=o 

+ E lk« +1 ) - ^(ui)||^ ( a) + E ll^u* 1 ) - d t (u{) 

1=0 1=0 

is uniformly bounded. For the first sum, we used the inequality 

ll m /i m h.llL2(f7) S « H v /illL2(f2)' 

see e.g. [1, 18]. We thus get 

N ~ 1 f tl+1 t — tt 

hhk - lhk\\h(n T ) = E / Ht£ + ~~ k~^ h+1 ~ 7 ^ ~ ^"^(n) 

<*£ll^WX (n) ^0, 

£=0 

13 



Il 2 (Q) 



and analogously 



i=0 Jtl 
N-l 



— Kfh ~ IhJ ~ Ih IIL2(Q) 



£=0 



We thus conclude that the limits lim^/^ = lim7^ fc coincide in L 2 (f2y). From the continuous 
inclusions H 1 (r2^) <e L 2 (H 1 ) C L 2 (fir) and the uniqueness of weak limits, we even conclude 
the convergence properties in L 2 (H 1 ), resp. in H 1 (f2j'). By use of the Weyl theorem, we may 
extract yet another subsequence to see pointwise convergence, i.e. (29d). From 

UN - lHi^nr) < IIH - |m^||| L 2 (nr) + || jm.- fe j - l|| L a(n r ), and 
lll m ftfe(^ -)l ~ 1 IIl 2 (^) < ^max||Vm^|| L 2 (n) 0, 

we finally conclude |m| = 1 almost everywhere in Qt, which is the desired result. □ 

In the remainder of this chapter, we will prove that the limiting tupel (m, u) is indeed a 
weak solution in the sense of Definition 1. First, we identify the limit function v with the 
time derivative of m. The following result can be found e.g. in [1]. 

Lemma 11. The limit function v £ L 2 (Qt) equals the time derivative of m, i.e. v = c^m 
almost everywhere in Oy. □ 

We have now collected all ingredients for the proof of our main theorem. 

Proof of Theorem 4. Let (C,i/>) £ C°°(fi T ) x C c °° ([0, T); C°°(fi)) be arbitrary. We define 
testfunctions by ((p h ,if) h )(t, ■) := (Z h (m.y ik x £),I h iJ), ) (t, ■). With the notation from above, 
we integrate equation (16) in time to obtain 

« / ( v ^> <Ph) + / ((™hk x v m)> <Ph) = ~°e I (V(m^ + ekvu), V^)) 
Jo Jo Jo 

+ / (h m (u^,m^),y> fc ) - / (7r(m^),y> h ). 
Jo Jo 

Then the magnetostrictive component is again given by 

[Mu^,m- )]< := \™ pg a hk v(m- hk ) p , with cr" = \ e (e(u~ hk ) - e m {m^)). 
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The definition (p h (t, •) := Ih(m hk x £)(£, •) and the approximation properties of the nodal 
interpolation operator show 

/ ((« v m + m hk x ^hk)X m hk x 0) + k 9 [ (Vv~ fc , V(m^ x C)) 
Jo Jo 

+ C e f (Vm tt) V(mJi x C)) 
jo 

- / (h m (u~ fe , m- fe ), (m- fe x Q) 
Jo 

+ / (^(m^),(m^ x C)) 
Jo 
= 0(h). 

Passing to the limit and using the strong L 2 (nr)-convergence of (m^ fc x £) towards (m x £), 
we get 

/ (i av hk + m hk x v hk)^ m hk x 0) — >[ (K + mx m t ),(mx C)), 
jo Jo 

kO [ (Vv^, V(m^ x C)) — ► 0, and 
Jo 

/ (Vm u ,VKx())^/ (Vm,V(mxC)), 
Jo JO 

as (h,k) (0,0), cf. [1, Proof of Thm. 2]. Here, we have used the boundedness of 

^ll^ v hJI:L 2 (n t ) wn i cn follows from 9 G (1/2, 1], see Corollary 9. Next, the weak convergence 

of 7r(m^ fe ) from (21) yields 

/ (t(%),Kx 0) —> [ (7r(m),(mx 0)- 

JO JO 

As for the magnetostrictive component, we have to show 

h m (% fc , m- fc ) h m (u, m) in L 2 (fi T ), 

where it obviously suffices to show the desired property componentwise. With the definition 
from (6) and the boundedness of A m , a direct computation proves 

||s m (m,- fc ) - e-(m)|| 2 L2(nT) < \\m hk - m \\l H n T ), 

cf. [24]. The pointwise convergence of m^ fc from (29d) in combination with Lebesgue's 
dominated convergence theorem, now yield the strong convergence m^, • £ — > m • £. For any 
indices i, j,p = 1, 2, 3 this shows 

((£ m (™hk))v( m hk)pX P ) ->• (e m (m) lJ m P X P )- 
By definition of the magnetostrictive component it therefore only remains to show 

( e ( u hk))v( m hk)p ~* £ ( u )ij ■ m P in L 2 (^r) 
for any combination i,j,p = 1, 2, 3 of indices. Analogously to above, this can be seen by 

(( £ ( u hk))ij( m hk) P X P ) = ( £ (^hk)ij,(^hk)pC P ) -> (^(u) lJ ,m p C p ) = (£(u) 4j m p , C p ) 

15 



for all £ G C°°(Qt), where the convergence of £{u hk ) towards e(u) particularly follows from 
the convergence of u^ k towards u in L 2 (H 1 ). So far, we have thus proved 

/ ((am t + m x m t ), (m x £)) = — C e / (Vm, V(m x £)) 
Jo Jo 

+ [ (h m (u,m),(m x C)) - / (ff(m),(mxC)). 
Jo Jo 

Proceeding as in [5], we conclude (11) by use of elementary pointwise calculations. The 

equality m(0, ■) = m° in the trace sense follows from the weak convergence — ^ m 

in H 1 (fix') and thus weak convergence of the traces. The equality u(0, ■) = u° follows 

analogously. 

To prove (12), we argue similarly. From (17), we obtain 

[ T ((u hk ) t ,i> h )+ / T (A e 6(u+),6(^))= / T (A e £ m (m+),6(^)). 
Jo Jo Jo 

For the first summand on the left-hand side, we perform integration by parts in time and, 

due to the shape of if) h , get 

f ((u hk ) t ,ip h ) = - [ (u hfc ,(^*) + Wr,o,^ h (r,o)-(uhfc(o, 0,^(0, •))• 

JO JO V * ' ' v ' 



Passing to the limit (h, k) — > 0, we see 



/ {(u hk ) u il> h ) — [ (u,^)-(u(0,-),V(0,-))- 
Jo Jo 



Here, we have used the assumed convergence of the initial data. From (29h) and the definition 
of u^" fc , we get u^ k = dtUhk, and therefore, by use of weak lower semi- continuity, conclude 

l|u - ^ u IIl2(h) < hminf ||u+ - d t u hk \\l 2{n) = 
whence u = d t u almost everywhere in fiy The convergence of the terms 

/ (A e £ (u+),e(^)) — ► / T (A e e( U ), £(</>)) and 
Jo Jo 

/ T (A e £ -(m+),e(^)) — ► / T (A e £ m (m),6W) 
Jo Jo 

is straightforward. In summary, we have thus shown (12). 

It remains to show the energy estimate (13). From the discrete energy estimates (25) 

and (27), in combination with Korn's inequality, we get for any t' G [0, T] with t' G [te,te+i) 

!|Vm+ (011^(0) + IKJ^,) + ||Vu+ (tOll^cn) + \\^tk^)\\b(n) 

= ||Vm+ (Oll^n) + / \K k (t)\\h(n) + HVu+ (t')\\l {n) + ||u+ (f 

Jo 



< ||Vm+ (Oll^nj + / IK^IIl^) + l|Vu+ (0 \\h m + \K k (t')\\h ( 



<c, 
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for some constant C which is independent of h and k. Integration in time thus yields for any 
measurable set TC [0, T] 

J II Vm+ (011^(0) + J hZAfto) + J l|Vu+ (Oilmen) + J 11^(011^(0) < jc. 
Weak semi-continuity now shows 

Vm|£ 2(n) + / ||m t ||^ 2( o ;) + / ||Vu||£ 3(n) + / ||u t ||^ (n) < C, 
% J% Ji J% J% 

which concludes the proof. □ 

Remark. Finally, we like to comment on the choice of < 9 < 1/2 to emphasize how 
variants of Theorem 4 are read and proved. 

(1) For < 9 < 1/2, one has to bound the term k 2 Yle=o II^ v 1IIl2(q) ^ n Corollary 9 in 
order to prove boundedness of the discrete quantities. This can be achieved by using 
an inverse estimate HVv^H^q) < 7^II v /iIIl2(q)- ^ e ^ a ^ er term can then be absorbed 
into the term fc^^~ 1 1 1 1 l 2 (s^) an< ^ ^ us convergence as p- — > 0. 

(2) For the intermediate case 9 — | ; the limit k9 (Vv^, V(m^ fc x ^)) — y is no longer 
valid, see [1, Proof of Thm. 2]. As suggested in [1], this can be circumvented by an 
inverse estimate provided that | — > 0. 

□ 



6. Numerical experiments 

In this section, we investigate the performance of the proposed Algorithm 2 empirically. 
We compare Algorithm 2 with the midpoint scheme from [9, 25] for the following blow-up 
benchmark example in 2D, proposed in [7, 8]: On Q = (— 0.5,0. 5) 2 , we solve problem (9). 
For some fixed parameter s6l and A := (1 — 2|x|) 4 /s , the initial magnetization reads 

f (0,0,-1) for |x| > 0.5 

mo(x> = i ^mp w ^ ws. 

This choice of mo is motivated in [8] in order to form some singularity at the center x = 0. 

The triangulation % used in the numerical simulation is defined through a positive integer 
r and consists of 2 2r+1 halved squares with edge length h = 2~ r . Since this triangulation is of 
Delaunay type, the angle condition (20) is satisfied, see [1] . For the magnetostriction, we took 
A e and A m to be 2 x 2 tensors with A e = Af m = X e 2222 = C e , X m = Xf in = A^ 222 = C m . 
for given constants C e , C m > 0. For comparison, the midpoint scheme is computed as 
described in [25]. The latter requires a fixed point iteration, where the tolerance is chosen as 
e = 10 -10 , and a time-step k m > that satisfies the associated mesh size condition k m < Ch 2 
from [4, 9, 25]. In our computations for the midpoint scheme, we thus chose k m = h 2 /10. The 
time-step size of Algorithm 2 is denoted by k a and does not depend on the spatial mesh-size, 
i.e. as h is decreased, k a remains fixed throughout the computations while k m needs to be 
stepwise decreased. The linear systems of Algorithm 2 are solved using an iterative method, 
and the constraint on the space JC m e is incorporated via the Lagrange multiplier approach 
from [17, 18]. 
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14157 
N.A.N 


0.059 
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Table 1. Comparison of CPU and blow-up time Tb for Algorithm 2 and 
midpoint scheme [9] with a — 1, T=0.3[s], and C e = C m = 0. 




0.01 0.02 0.03 0.04 0.05 0.06 



Time[s] 

Figure 1. W l, °°(Q,) semi- norm from Algorithm 2 for different values of h, 
with C e = C m = r=5, s=4, k a = 10"5. 

In a first experiment, we consider the exchange-only case of LLG and thus neglect magne- 
tostrictive effects as well as all other field contributions, i.e. 7r(-) = 0. We compare the two 
algorithms for a = 1, s = 4, 9 = 1, k a = 10" 5 , C e = C m = 0, and T = 0.3[s]. As initital value 
for k m , we choose 10 -5 . In Table 1, we investigate the overall computation time as well as 
the empirical blow-up time, which is the time when the |m(i)|i i00 = ||Vm(i)||£°o(n) reaches 
its maximum. We observe that Algorithm 2 leads to significantly lower computation time 
than the midpoint scheme. Moreover, the blow-up time seems to vary between the different 
schemes. We observed that the blow-up times can be brought more in line with each other 
if the time-step size k a is a little decreased (not displayed). As expected, in comparison to 
the midpoint scheme, Algorithm 2 can work with larger time-steps for reasonably small h. 

On the other hand, the blow-up time seems to increase as h becomes smaller as shown in 
Figure 1. This effect was not observable in [9] due to the fixed-point iteration and thus the 
coupling of h and k m . Our empirical observation raises the question of the mere existence 
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Algorithm 4.1 (alpha=1/1000) 




Time[s] 



Figure 2. Evolution of the Energy for different values of a, with C e = C m = 
r=5, s=4, k a = 10-5, k m = 10~5. 

and behaviour of the blow-up time in case of exchange only. Put explicitly, the finite time 
blow-up might be a numerical artifact stemming from insufficient spatial resolution. 

In Figure 2, we plot the discrete energy for the exchange-only case with C e = C m = de- 
fined by E(m, t) = !||Vm(i)||£ 2 (m> anc ^ investigate the stability of the respective algorithms 
when a becomes small. The kinks in the graph coincide with the empirical blow-up times. 
We observe that both algorithms seem to be stable, while Algorithm 2 provides a stronger 
energy decay for small values of a. 
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2899 
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0.0074 
N.A.N 



Table 2. Comparison of CPU and blow-up time Tg for Algorithm 2 and 
midpoint scheme [25] with a = 1, T=0.015[s], C e = 40 and C m = 10. 



In a second experiment, we include magnetostriction. In Table 2, we compare the two 
algorithms for C e = 40, C m = 10, and T = 0.015[s] and observe similar results as for the 
exchange-only case for LLC Even with included magnetostrictive effects, the blow-up time 
varies a lot depending on the spatial resolution h, as well as between the different schemes. 
Finally, in the Figures 3 and 4, we compare the computational results of the two different 
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Figure 3. Evolution of HmJ^j = 1,3 with C e = C m = a = 1/64, r=5, 
s=4, k a = k m = 10"5. 

schemes. Due to the symmetry of the problem, we only show the evolution of the L 2 (f2)- 
average 

mj\\ L 2 {n) , j = 1,3 

for the mi and components of the magnetization. 

Overall, we conclude that the results of our algorithm are in good agreement with the 
midpoint scheme, more feasible for small a, throughout much faster to compute, and finally 
easier to implement. 
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